En la sesión de hoy veremos:

Regresión lineal simple

Los siguientes son datos del libro: “Understandable Statistics”. Los datos son para seguros de autos. Cada renglón es una región geográfica de Suecia.

library(readxl)
datos <- read_xls("../R_docentes/datos/slr06.xls")
head(datos)
## # A tibble: 6 × 2
##       X     Y
##   <dbl> <dbl>
## 1   108 392. 
## 2    19  46.2
## 3    13  15.7
## 4   124 422. 
## 5    40 119. 
## 6    57 171.
dim(datos)  # dimensión del dataframe
## [1] 63  2
colnames(datos) <- c("n_claims","total_payment")

Estadísticas descriptivas

summary(datos)
##     n_claims     total_payment   
##  Min.   :  0.0   Min.   :  0.00  
##  1st Qu.:  7.5   1st Qu.: 38.85  
##  Median : 14.0   Median : 73.40  
##  Mean   : 22.9   Mean   : 98.19  
##  3rd Qu.: 29.0   3rd Qu.:140.00  
##  Max.   :124.0   Max.   :422.20
cor(datos)
##                n_claims total_payment
## n_claims      1.0000000     0.9128782
## total_payment 0.9128782     1.0000000
# Podemos generar una gráfica de la matriz de correlaciones con el paquete corrplot
library(corrplot)
## corrplot 0.92 loaded
corrplot(cor(datos), method = "ellipse")

# gráfica de dispersión de puntos
plot(datos$n_claims, datos$total_payment)
abline(lm(total_payment ~ n_claims, data = datos)) # Agrega la línea de mínimos cuadrados
# regresión no paramétrica
a <- loess(total_payment ~ n_claims, data = datos) # Regresión no paramétrica
j <- order(a$x)
lines(a$x[j], a$fitted[j], col = "red", lwd=3)

# Versión ggplot:
library(ggplot2)

ggplot(datos,aes(n_claims, total_payment)) +
   geom_point(alpha = 0.5) +
   geom_smooth(method = "lm", se = F) +
   geom_smooth(method = "loess", se = F, color = "red")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Ajuste del modelo de regresión

Ya lo hicimos arriba, para la gráfica, pero vamos a ver los detalles.

mod1 <- lm(total_payment ~ n_claims, data = datos)
mod1 # devuelve los parámetros estimados
## 
## Call:
## lm(formula = total_payment ~ n_claims, data = datos)
## 
## Coefficients:
## (Intercept)     n_claims  
##      19.994        3.414
summary(mod1)  # nos da todos los detalles: estimación, errores estándar, anova. 
## 
## Call:
## lm(formula = total_payment ~ n_claims, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.561 -24.051  -0.347  23.432  83.977 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.9945     6.3678    3.14   0.0026 ** 
## n_claims      3.4138     0.1955   17.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.94 on 61 degrees of freedom
## Multiple R-squared:  0.8333, Adjusted R-squared:  0.8306 
## F-statistic:   305 on 1 and 61 DF,  p-value: < 2.2e-16
model.matrix(mod1) # muestra la matriz de diseño X correspondiente al modelo
##    (Intercept) n_claims
## 1            1      108
## 2            1       19
## 3            1       13
## 4            1      124
## 5            1       40
## 6            1       57
## 7            1       23
## 8            1       14
## 9            1       45
## 10           1       10
## 11           1        5
## 12           1       48
## 13           1       11
## 14           1       23
## 15           1        7
## 16           1        2
## 17           1       24
## 18           1        6
## 19           1        3
## 20           1       23
## 21           1        6
## 22           1        9
## 23           1        9
## 24           1        3
## 25           1       29
## 26           1        7
## 27           1        4
## 28           1       20
## 29           1        7
## 30           1        4
## 31           1        0
## 32           1       25
## 33           1        6
## 34           1        5
## 35           1       22
## 36           1       11
## 37           1       61
## 38           1       12
## 39           1        4
## 40           1       16
## 41           1       13
## 42           1       60
## 43           1       41
## 44           1       37
## 45           1       55
## 46           1       41
## 47           1       11
## 48           1       27
## 49           1        8
## 50           1        3
## 51           1       17
## 52           1       13
## 53           1       13
## 54           1       15
## 55           1        8
## 56           1       29
## 57           1       30
## 58           1       24
## 59           1        9
## 60           1       31
## 61           1       14
## 62           1       53
## 63           1       26
## attr(,"assign")
## [1] 0 1

Datos categóricos:

library(ggplot2)
datos2 <- iris[,c(1,5)] # Nos quedamos con las columnas uno y cinco
ggplot(datos2, aes(Sepal.Length)) +
       geom_histogram(bins = 10) + 
       facet_wrap(vars(Species))

Estadísticas por categoría:

# tradicional
tapply(datos2$Sepal.Length,datos2$Species, summary)
## $setosa
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.300   4.800   5.000   5.006   5.200   5.800 
## 
## $versicolor
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.900   5.600   5.900   5.936   6.300   7.000 
## 
## $virginica
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.900   6.225   6.500   6.588   6.900   7.900
# tidy
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
datos2 %>% 
       group_by(Species) %>%
       summarize(media = mean(Sepal.Length),
                 sd = sd(Sepal.Length),
                 n = n())
## # A tibble: 3 × 4
##   Species    media    sd     n
##   <fct>      <dbl> <dbl> <int>
## 1 setosa      5.01 0.352    50
## 2 versicolor  5.94 0.516    50
## 3 virginica   6.59 0.636    50

El modelo de regresión en este caso:

# Aquí R está considerando a Species como un factor, aunque es sólo un vector de caracteres. 
# Se recomienda convertir Species a factor en caso de problemas. 
mod2 <- lm(Sepal.Length ~ Species, data = datos2)
summary(mod2)
## 
## Call:
## lm(formula = Sepal.Length ~ Species, data = datos2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6880 -0.3285 -0.0060  0.3120  1.3120 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.0060     0.0728  68.762  < 2e-16 ***
## Speciesversicolor   0.9300     0.1030   9.033 8.77e-16 ***
## Speciesvirginica    1.5820     0.1030  15.366  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared:  0.6187, Adjusted R-squared:  0.6135 
## F-statistic: 119.3 on 2 and 147 DF,  p-value: < 2.2e-16
X <- model.matrix(mod2)
head(X)
##   (Intercept) Speciesversicolor Speciesvirginica
## 1           1                 0                0
## 2           1                 0                0
## 3           1                 0                0
## 4           1                 0                0
## 5           1                 0                0
## 6           1                 0                0
# También se puede quitar la ordenada del modelo
mod2a <- lm(Sepal.Length ~ Species + 0, data = datos2)
mod2a
## 
## Call:
## lm(formula = Sepal.Length ~ Species + 0, data = datos2)
## 
## Coefficients:
##     Speciessetosa  Speciesversicolor   Speciesvirginica  
##             5.006              5.936              6.588

Para extraer información del objeto

coefficients(mod1)  # extrae los coeficientes
## (Intercept)    n_claims 
##   19.994486    3.413824
fitted(mod1)        # extrae los valores ajustados
##         1         2         3         4         5         6         7         8 
## 388.68743  84.85713  64.37419 443.30861 156.54743 214.58243  98.51243  67.78802 
##         9        10        11        12        13        14        15        16 
## 173.61655  54.13272  37.06360 183.85802  57.54654  98.51243  43.89125  26.82213 
##        17        18        19        20        21        22        23        24 
## 101.92625  40.47743  30.23596  98.51243  40.47743  50.71890  50.71890  30.23596 
##        25        26        27        28        29        30        31        32 
## 118.99537  43.89125  33.64978  88.27096  43.89125  33.64978  19.99449 105.34007 
##        33        34        35        36        37        38        39        40 
##  40.47743  37.06360  95.09860  57.54654 228.23772  60.96037  33.64978  74.61566 
##        41        42        43        44        45        46        47        48 
##  64.37419 224.82390 159.96125 146.30596 207.75478 159.96125  57.54654 112.16772 
##        49        50        51        52        53        54        55        56 
##  47.30507  30.23596  78.02949  64.37419  64.37419  71.20184  47.30507 118.99537 
##        57        58        59        60        61        62        63 
## 122.40919 101.92625  50.71890 125.82302  67.78802 200.92713 108.75390
residuals(mod1)     # extrae los residuales del modelo (varios tipos de residual, con opción type = )
##           1           2           3           4           5           6 
##   3.8125698 -38.6571334 -48.6741920 -21.1086072 -37.1474282 -43.6824287 
##           7           8           9          10          11          12 
## -41.6124276   9.7119844  40.3834540  11.1672786 -16.1636036  64.2419834 
##          13          14          15          16          17          18 
## -34.0465449 -58.9124276   4.9087493 -20.2221329  32.9737488  10.4225729 
##          19          20          21          22          23          24 
## -25.8359564  14.4875724 -25.6774271  -2.0188978   1.3811022 -17.0359564 
##          25          26          27          28          29          30 
## -15.0953690  33.6087493 -21.8497800   9.8290430 -15.9912507   4.4502200 
##          31          32          33          34          35          36 
## -19.9944858 -36.1400748 -25.8774271   3.2363964  66.4013959  -0.3465449 
##          37          38          39          40          41          42 
## -10.6377229  -2.8603685 -21.0497800 -15.0156627  25.5258080 -22.4238994 
##          43          44          45          46          47          48 
##  21.3387483   6.4940425 -44.9547816 -86.5612517 -36.2465449 -19.5677219 
##          49          50          51          52          53          54 
##  28.7949258   9.6640436  64.0705137  28.6258080 -32.4741920 -39.1018392 
##          55          56          57          58          59          60 
##   8.2949258  14.3046310  72.0908074  35.9737488  36.6811022  83.9769839 
##          61          62          63 
##  27.7119844  43.6728656  78.7461017
cooks.distance(mod1) # distancias de Cook para diagnósticos
##            1            2            3            4            5            6 
## 2.183354e-03 9.758274e-03 1.788235e-02 1.180318e-01 1.376184e-02 4.115139e-02 
##            7            8            9           10           11           12 
## 1.098506e-02 6.900603e-04 2.034982e-02 1.047074e-03 2.699222e-03 5.912148e-02 
##           13           14           15           16           17           18 
## 9.375163e-03 2.201761e-02 2.283703e-04 4.832807e-03 6.913315e-03 1.074464e-03 
##           19           20           21           22           23           24 
## 7.539119e-03 1.331520e-03 6.521462e-03 3.558483e-05 1.665282e-05 3.277967e-03 
##           25           26           27           28           29           30 
## 1.549092e-03 1.070539e-02 5.155686e-03 6.228258e-04 2.423610e-03 2.138727e-04 
##           31           32           33           34           35           36 
## 5.177366e-03 8.355639e-03 6.623448e-03 1.082146e-04 2.801473e-02 9.712964e-07 
##           37           38           39           40           41           42 
## 2.907237e-03 6.385946e-05 4.785060e-03 1.561876e-03 4.917973e-03 1.237114e-02 
##           43           44           45           46           47           48 
## 4.743961e-03 3.709949e-04 3.985862e-02 7.806390e-02 1.062591e-02 2.507445e-03 
##           49           50           51           52           53           54 
## 7.537632e-03 1.054847e-03 2.779156e-02 6.185042e-03 7.959831e-03 1.086948e-02 
##           55           56           57           58           59           60 
## 6.255010e-04 1.391051e-03 3.617134e-02 8.228508e-03 1.174684e-02 5.039858e-02 
##           61           62           63 
## 5.618318e-03 3.436581e-02 4.006284e-02
# Un poco más avanzado:
library(broom)
tidy(mod1) # regresa las estadísticas en forma de dataframe para programación
## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    20.0      6.37       3.14 2.60e- 3
## 2 n_claims        3.41     0.195     17.5  2.05e-25
augment(mod1) # un dataframe con varias estadísticas asociadas a cada par de datos. 
## # A tibble: 63 × 8
##    total_payment n_claims .fitted .resid   .hat .sigma  .cooksd .std.resid
##            <dbl>    <dbl>   <dbl>  <dbl>  <dbl>  <dbl>    <dbl>      <dbl>
##  1         392.       108   389.    3.81 0.230    36.2 0.00218       0.121
##  2          46.2       19    84.9 -38.7  0.0163   35.9 0.00976      -1.08 
##  3          15.7       13    64.4 -48.7  0.0188   35.7 0.0179       -1.37 
##  4         422.       124   443.  -21.1  0.318    36.1 0.118        -0.711
##  5         119.        40   157.  -37.1  0.0245   35.9 0.0138       -1.05 
##  6         171.        57   215.  -43.7  0.0503   35.8 0.0412       -1.25 
##  7          56.9       23    98.5 -41.6  0.0159   35.8 0.0110       -1.17 
##  8          77.5       14    67.8   9.71 0.0182   36.2 0.000690      0.273
##  9         214         45   174.   40.4  0.0303   35.9 0.0203        1.14 
## 10          65.3       10    54.1  11.2  0.0208   36.2 0.00105       0.314
## # … with 53 more rows
glance(mod1)  # estadísticas generales del modelo
## # A tibble: 1 × 12
##   r.squared adj.r.squa…¹ sigma stati…²  p.value    df logLik   AIC   BIC devia…³
##       <dbl>        <dbl> <dbl>   <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>
## 1     0.833        0.831  35.9    305. 2.05e-25     1  -314.  634.  640.  78797.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## #   variable names ¹​adj.r.squared, ²​statistic, ³​deviance

Predicciones

Usando el primer modelo con variables continuas. Necesitamos crear un arreglo con los valores de los predictores que queremos considerar en la evaluación

ncasos <- data.frame(n_claims=80:100)
predict(mod1,ncasos)
##        1        2        3        4        5        6        7        8 
## 293.1004 296.5142 299.9280 303.3418 306.7557 310.1695 313.5833 316.9971 
##        9       10       11       12       13       14       15       16 
## 320.4110 323.8248 327.2386 330.6524 334.0663 337.4801 340.8939 344.3077 
##       17       18       19       20       21 
## 347.7215 351.1354 354.5492 357.9630 361.3768
# Serìa mejor tener en un sólo arreglo los valores nuevos y sus predicciones
predicciones <- ncasos %>%
                mutate(total_payment = predict(mod1,ncasos))
# Podemos agregar las predicciones a una gráfica
plot(datos$n_claims, datos$total_payment)
abline(mod1)
points(predicciones, col = "red", pch = 19)

# Con ggplot
datos %>% 
  ggplot(aes(x = n_claims, y = total_payment)) +
  geom_point() + 
  geom_smooth(method = "lm", se = F) +
  geom_point(
      data = predicciones,
      color = "red"
  )
## `geom_smooth()` using formula = 'y ~ x'

También podemos obtener los intervalos de confianza para las predicciones. Hay dos opciones: una para la estimación de un valor promedio, y la otra para una nueva observación.

A <- as.data.frame(predict(mod1,interval = "confidence",level=0.8))
B <- as.data.frame(predict(mod1,interval = "prediction",level=0.8))
## Warning in predict.lm(mod1, interval = "prediction", level = 0.8): predictions on current data refer to _future_ responses
# añade intervalos de confianza a la gráfica
plot(datos$n_claims,datos$total_payment)
abline(mod1)
lines(datos$n_claims, A$lwr, col = "red")
lines(datos$n_claims, A$upr, col = "red")
lines(datos$n_claims, B$lwr, col =  "green")
lines(datos$n_claims, B$upr, col = "green")

# Con ggplot:
datos %>% 
  ggplot(aes(x = n_claims, y = total_payment)) +
  geom_point() + 
  geom_smooth(method = "lm", se = T) +
  geom_line(aes(y=B$lwr), color = "orange", linetype = "dashed") +
  geom_line(aes(y=B$upr), color = "orange", linetype = "dashed") 
## `geom_smooth()` using formula = 'y ~ x'

Diagnósticos para el modelo

Hay una serie de gráficas para verificar los supuestos del modelo

par(mfrow = c(2,2)) # crea un arreglo para múltiples gráficas
plot(mod1) # gráficas de diagnóstico

Transformación de variables

Consideremos los siguientes datos: es un registro de 7 especies comunes de peces en ventas de mercado. Los datos son de Kaggle. Seleccionamos un grupo particular, el grupo de lubinas (perch)

Fish <- read.csv("https://raw.githubusercontent.com/jvega68/EA3/master/datos/Fish.csv")
unique(Fish$Species)
## [1] "Bream"     "Roach"     "Whitefish" "Parkki"    "Perch"     "Pike"     
## [7] "Smelt"
perch <- Fish %>%
         filter(Species == "Perch")
with(perch, plot(Length1,Weight)) # la relación no es lineal 

with(perch, plot(Length1^3,Weight)) # Vemos que la relación se lineariza con esta transformación

Para hacer la transformación en la regresión, podemos aplicar funciones a los predictores, pero en el caso de potencias, requerimos añadir el operador I(), porque el símbolo de potencia ^ tiene una función especial, que se verá más adelante

(mod3 <- lm(Weight ~ I(Length1^3), data = perch))
## 
## Call:
## lm(formula = Weight ~ I(Length1^3), data = perch)
## 
## Coefficients:
##  (Intercept)  I(Length1^3)  
##      -0.1175        0.0168

Para realizar la predicción, no es necesario elevar los datos al cubo, eso se hace como parte del modelo:

nuevos_datos <- data.frame( Length1 = seq(10,40,5))
datos_predicc <- data.frame(x = nuevos_datos, 
                            y = predict(mod3, nuevos_datos))
with(perch, plot(Length1^3,Weight)) # Vemos que la relación se lineariza con esta transformación
abline(mod3)
# agrega los puntos rehaciendo la transformación para que los grafique correctamente. 
points(datos_predicc$Length1^3, datos_predicc$y, col = "blue", pch = 19, cex = 2)

El tratamiento para transformar la respuesta es un poco diferente. Los siguientes datos también son de Kaggle y tienen información de tres campañas de marketing. Cada renglón es un anuncio.

Nos fijamos en tres variables: - spent: gasto en USD - impressions: las veces que las personas vieron los anuncios - clicks: las veces que las personas dieron clcks en los anuncios

fb_ad <- read.csv("https://raw.githubusercontent.com/jvega68/EA3/master/datos/fb_adv_data.csv")
str(fb_ad)
## 'data.frame':    1143 obs. of  15 variables:
##  $ ad_id              : int  708746 708749 708771 708815 708818 708820 708889 708895 708953 708958 ...
##  $ reporting_start    : chr  "17/08/2017" "17/08/2017" "17/08/2017" "30/08/2017" ...
##  $ reporting_end      : chr  "17/08/2017" "17/08/2017" "17/08/2017" "30/08/2017" ...
##  $ campaign_id        : chr  "916" "916" "916" "916" ...
##  $ fb_campaign_id     : chr  "103916" "103917" "103920" "103928" ...
##  $ age                : chr  "30-34" "30-34" "30-34" "30-34" ...
##  $ gender             : chr  "M" "M" "M" "M" ...
##  $ interest1          : int  15 16 20 28 28 29 15 16 27 28 ...
##  $ interest2          : int  17 19 25 32 33 30 16 20 31 32 ...
##  $ interest3          : int  17 21 22 32 32 30 17 18 31 31 ...
##  $ impressions        : num  7350 17861 693 4259 4133 ...
##  $ clicks             : int  1 2 0 1 1 0 3 1 1 3 ...
##  $ spent              : num  1.43 1.82 0 1.25 1.29 ...
##  $ total_conversion   : int  2 2 1 1 1 1 1 1 1 1 ...
##  $ approved_conversion: int  1 0 0 0 1 1 0 1 0 0 ...
with(fb_ad, plot(spent, impressions)) # los puntos están muy concentrados en la parte baja de la gráfica

with(fb_ad, plot(sqrt(spent), sqrt(impressions))) # se aprecian mejor los datos de la parte baja

# filtramos los casos con 0 impresiones o con 0 gasto
fb1 <- fb_ad %>% 
         filter(!((spent < 10) | (impressions == 0)))
dim(fb1)
## [1] 263  15
with(fb1, plot(sqrt(spent), sqrt(impressions))) # se aprecian mejor los datos de la parte baja

mod4 <- lm(sqrt(impressions) ~ sqrt(spent), data = fb1)
mod4
## 
## Call:
## lm(formula = sqrt(impressions) ~ sqrt(spent), data = fb1)
## 
## Coefficients:
## (Intercept)  sqrt(spent)  
##      -32.06        65.98
# prediccion
nvos_datos <- data.frame(spent = seq(0, 600, 100))
# la predicción devuelve los valores de la variable predictora en raíz cuadrada, por lo que necesitamos regresar
# a la escala original
datos_pred <- data.frame(
              spent = seq(0, 600, 100),
              sqrt_impressions = predict(mod4, nvos_datos))
datos_pred$impressions = datos_pred$sqrt_impressions^2 # agrega los datos en la escala original
plot(fb1$spent, fb1$impressions)
points(datos_pred$spent, datos_pred$impressions, col = "orange",pch = 19, cex = 2)

Regresión lineal múltiple

Revisaremos algunos puntos importantes en el análisis de regresión lineal múltiple, particularmente las diferencias de RLM con RLS. Muchas de las cosas que revisamos en RLS se aplican directamente en RLM.

Especificación de modelos en R en RLM

A continuación se muestra notación de modelos multivariados.

  • y ~ x \(y=b_0+b_1x\) Modelo de regresión lineal simple
  • y ~ -1 + x \(y= b_1x\) Modelo sin ordenada al origen
  • y ~ x + I(x^2) \(y=b_0 + b_1x+b_2x^2\) Modelo polinomial.
  • y ~ x1 + x2 \(y = b_0 + b_1x_1 + b_2x_2\) Modelo de primer orden sin interacción.
  • y ~ x1:x2 \(y = b_0 + b_1x_1x_2\) Modelo de interacción de primer orden
  • y ~ x1*x2 \(y = b_0 + b_1x_1 + b_2x_2 + b_3x_1x_2\) Modelo de primer orden completo. Un modelo equivalente es: - y ~ x1 + x2 + x1:x2
  • y ~ (x1 + x2 + x3)^2 \(y = b_0 + b_1x_1 + b_2x_2 + b_3x_3 +b_4x_1x_2 +b_5x_2x_3 +b_6x_1x_3\) Interacciones dobles

Consideremos el siguiente ejemplo con el archivo de Fish

Bream <- subset(Fish, Species == "Bream", select = -Species) # otra forma de tomar subconjuntos de datos.
mod5 <- lm(Weight ~ Length3 + Height + Width, data = Bream)
summary(mod5)
## 
## Call:
## lm(formula = Weight ~ Length3 + Height + Width, data = Bream)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -169.057  -26.519   -1.592   29.902  109.631 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1046.601     90.341 -11.585 8.54e-13 ***
## Length3        12.033      7.615   1.580 0.124248    
## Height         62.820     16.685   3.765 0.000699 ***
## Width          45.898     35.415   1.296 0.204548    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 52.76 on 31 degrees of freedom
## Multiple R-squared:  0.942,  Adjusted R-squared:  0.9364 
## F-statistic: 167.8 on 3 and 31 DF,  p-value: < 2.2e-16
# Podemos aplicar directamente a los datos originales
lm(Weight ~ Length3 + Height + Width, data = Fish, subset = Species == "Bream")
## 
## Call:
## lm(formula = Weight ~ Length3 + Height + Width, data = Fish, 
##     subset = Species == "Bream")
## 
## Coefficients:
## (Intercept)      Length3       Height        Width  
##    -1046.60        12.03        62.82        45.90
lm(Weight ~ Length3 + Height + Width, data = Fish, subset = Species == "Roach")
## 
## Call:
## lm(formula = Weight ~ Length3 + Height + Width, data = Fish, 
##     subset = Species == "Roach")
## 
## Coefficients:
## (Intercept)      Length3       Height        Width  
##    -319.924        8.608       -1.875       73.701

Scatterplot matrix: datos en pares

plot(Fish[,-1], col = factor(Fish$Species))

# Equivalente en ggplot 
library(GGally)  # Perdón, omití decirles que instalaran este paquete
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
Fish %>% 
  group_by(Species) %>% 
  ggpairs(aes(color = Species), alpha = 0.4)
## Warning in warn_if_args_exist(list(...)): Extra arguments: 'alpha' are being
## ignored. If these are meant to be aesthetics, submit them using the 'mapping'
## variable within ggpairs with ggplot2::aes or ggplot2::aes_string.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Otro paquete recién descubierto
library(scatterPlotMatrix)
Fish$Species <- factor(Fish$Species)
scatterPlotMatrix(Fish,zAxisDim = "Species")

Gráfica de correlaciones:

corrplot(cor(Fish[,-1]), method = "ellipse", order = "hclust")

Para realizar selección de modelos de manera automática

library(MASS) # No requiere instalación, es de los paquetes básicos que se instalan junto con R
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
mod5 <- lm(Weight ~ (Length3 + Height + Width)^3, data = Bream)
mod5
## 
## Call:
## lm(formula = Weight ~ (Length3 + Height + Width)^3, data = Bream)
## 
## Coefficients:
##          (Intercept)               Length3                Height  
##             2593.952               -72.945              -384.968  
##                Width        Length3:Height         Length3:Width  
##             -121.926                11.110                 2.541  
##         Height:Width  Length3:Height:Width  
##               47.000                -1.122
mod6 <- stepAIC(mod5,                      # modelo inicial
                scope = list(lower = ~1),  # se puede dar un modelo upper y un lower
                direction = "both" )       # foward o backward
## Start:  AIC=287.34
## Weight ~ (Length3 + Height + Width)^3
## 
##                        Df Sum of Sq   RSS    AIC
## - Length3:Height:Width  1    1517.9 82985 285.99
## <none>                              81467 287.34
## 
## Step:  AIC=285.99
## Weight ~ Length3 + Height + Width + Length3:Height + Length3:Width + 
##     Height:Width
## 
##                  Df Sum of Sq   RSS    AIC
## - Height:Width    1      4.74 82990 283.99
## - Length3:Width   1   2103.13 85088 284.86
## - Length3:Height  1   2864.02 85849 285.18
## <none>                        82985 285.99
## 
## Step:  AIC=283.99
## Weight ~ Length3 + Height + Width + Length3:Height + Length3:Width
## 
##                  Df Sum of Sq   RSS    AIC
## - Length3:Width   1    2610.2 85600 283.07
## - Length3:Height  1    3145.5 86135 283.29
## <none>                        82990 283.99
## 
## Step:  AIC=283.07
## Weight ~ Length3 + Height + Width + Length3:Height
## 
##                  Df Sum of Sq   RSS    AIC
## - Length3:Height  1     705.7 86306 281.36
## - Width           1    4146.3 89746 282.73
## <none>                        85600 283.07
## 
## Step:  AIC=281.36
## Weight ~ Length3 + Height + Width
## 
##           Df Sum of Sq    RSS    AIC
## - Width    1      4676  90982 281.21
## <none>                  86306 281.36
## - Length3  1      6951  93256 282.07
## - Height   1     39467 125772 292.54
## 
## Step:  AIC=281.21
## Weight ~ Length3 + Height
## 
##           Df Sum of Sq    RSS    AIC
## <none>                  90982 281.21
## - Length3  1     12717 103699 283.79
## - Height   1     62192 153173 297.44

Comparación de modelos (siempre que estén anidados)

anova(mod5, mod6)
## Analysis of Variance Table
## 
## Model 1: Weight ~ (Length3 + Height + Width)^3
## Model 2: Weight ~ Length3 + Height
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     27 81467                           
## 2     32 90982 -5   -9514.6 0.6307 0.6779

Transformación de Box-Cox: Determina la transformación optima de la variable de respuesta a normalidad.

boxcox(mod5)

Referencias

Buenas referencias para modelos lineales hay muchas, dentro de ellas:

  • John Fox. Applied Regression Analysis and Generalized Linear Models, 3rd. Ed. Sage, 2016.
  • John Fox & Sandford Weisberg. An R Companion to Applied Regression, 2nd. Ed. Sage, 2011. Este libro tiene su propio paquete: car (companion to applied regression)
  • W. N. Venables & B. D. Ripley. Modern Applied Statistics with S, 4th. Ed. Springer, 2002. Este libro es la base del paquete MASS